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AN INITIAL INVESTIGATION INTO METHODS OF COMPUTING 


TRANSONIC AERODYNAMIC SENSITIVITY COEFFICIENTS 
I . Introduct i on 

This report covers the period from January 1, 1988 thru June 30, 
1988. The primary tasks during this were were to complete and formally 
report on the methods developed for computing aerodynamic sensitivity 
coefficients using the quas t -anal yt i cal approach. 

1 1 . Per sonne 1 

The staff associated with this project during the present reporting 
period were Dr. Leland A. Carlson, Principal Investigator, and Hesham El 
Banna, Graduate Research Assistant. 

111. Research Prooress 

As previously reported <Ref. 1), this initial effort has 
concentrated on developing the quas i -anal yt i cal approach for 
two-dimensional transonic flow. To keep the problem computationally 
efficient and straightforward, this initial study has only considered 
two-dimensional flow and has modeled the problem using the transonic 
small perturbation equation equation. During this reporting period, 
this initial development has been essentially completed and formally 
reported in the Master of Science thesis of Hesham El Banna. Mr. El 
Banna received his M.S. degree in May 1988, and is continuing his 
studies at Texas A&M towards a Ph.D. 
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Since Mr. El Banna's, thesis is an excellent summary of much of the 


work to date, it is included as an appendix o-f this report. In 
addition, a shortened version o-f this thesis has been submitted -for 
possible presentation at the 1989 AIAA Aerospace Sciences Meeting under 
applied aerodynamics. 


1U . Future Efforts 

As you are aware, the present project has been granted a no-cost 
extension until 31 December 1988. Currently studies are underway to 
compute aerodynamic sensitivity coe-f -f i c i en ts over a a series of 
transonic freestream Mach numbers in order to determine how rapidly the 
coefficients vary with Mach number and whether or not the present method 
detects or predicts these variations properly. In addition, new 
solution schemes for the quas i -anal yt i cal equation are in progress of 
being developed and tested. The objective, of course, is to decrease 
the computational time required to solve the quasi -anal yt i cal equation. 
Particular attention is currently being devoted to tri-diagonal 
iterative schemes which may carry over to three dimensions. These 
efforts will be reported in the next progress report. 

Finally, it is anticipated that a proposal will be submjitted to 
continue this work and to extend it to three dimensional transonic flows 
about wings. 
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ABSTRACT 

Numerical Computation of Aerodynamic Sensitivity 
Coefficients in the Transonic and Supersonic 
Regimes. (May 1988) 

Hesham M. Elbanna, B.S., Cairo University, Egypt 
Chair of Advisory Committee: Dr. Leland Carlson 

The quas i - analytical approach is developed to compute 
the aerodynamic sensitivity coefficients in the transonic 
and supersonic flight regimes. Initial investigation 
verifies the feasibility of this approach as applied to 
the transonic small perturbation residual expression. 
Results are compared to those obtained by the direct 
(finite difference) approach and both methods are 
evaluated to determine their computational efficiencies 
A Gauss-Seidel procedure is used to solve the large set of 
equations associated with the quas i - analytical approach 
On a medium grid, the quas i - analytical method is more 
efficient than the finite difference approach. However, on 
a fine grid, time comparisons are not as competitive. 
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INTRODUCTION 

Over the past few years, computational fluid dynamics 
has evolved rapidly as a result of the immense 
advancements in the computational field and the impact of 
the use of computers on obtaining numerical solutions to 
complex problems. Accordingly, researchers are now capable 
of calculating aerodynamic forces on wing - body - nac e 1 1 e - 
empennage configurations subject to subsonic or transonic 
flows. A next logical step would be to compute the 
sensitivity of these forces to configuration geometry 

In the transonic regime, one of the main difficulties 
facing the aircraft designer is the prediction of the 
aerodynamic loads. The difficulty is caused by the fact 
that in this regime even the most primitive representation 
of the aerodynamics must be described by a nonlinear 
equation or a set of equations. In addition, aerodynamic 
prediction in the transonic regime is extremely important 
since it is in this speed range that most civil aircraft 
maneuver. Consequently, the transonic regime is probably 
the most critical flow regime for present day aircraft 


Format in accordance with AIAA J ournal . 



In order to improve the design of transonic vehicles, 
design codes are being developed which use optimization 
techniques, and, in order to be successful, these codes 
require aerodynamic sensitivity coefficients, which are 
defined as the derivatives of the aerodynamic functions 
with respect to the design variables. Obviously, it is 
desirable that such sensitivity coefficients be easily 
ob t a ine d . 

Consequently, the primary objective of this effort is 
to investigate the feasibility of using the quasi- 
analytical method^*^ for calculating the aerodynamic 
sensitivity derivatives in the transonic and supersonic 
flight regimes. As part of this work, the resulting 
sensitivity coefficients are compared to those obtained 
from the finite difference approach. Finally, both methods 
are evaluated to determine their computational 
efficiencies . 

As mentioned earlier, knowledge of the sensitivity 
coefficients is essential information in any design 
optimization process. Obviously these calculations cannot 
be performed without the availability of solutions for the 
problem under consideration. 

In the transonic regime, a variety of methods for 
computing solutions to the flow field do exist. These 



range from full Navier - S tokes solvers to transonic small 
perturbation equation solvers. The complexity of the 
equations that need to be solved depends upon the flow 
phenomena in question and the objective of the analysis 
Since it is not the objective of this work to develop 
flowfield algorithms, the present research uses the 
transonic small perturbation equation to determine and, 
mentioned earlier, verify the existence of efficient 
methods for calculating the aerodynamic sensitivity 
derivatives . 

This research is original in that it aims to include 
sensitivity analysis procedures as part of aerodynamic 
analyses. Thus, it will provide a reference point for 
aeronautical engineers^ who need sensitivity information 
when conducting aerodynamic optimization as part of the 
aircraft design process. 
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BACKGROUND 

Most recently, sensitivity methodology has been 
successfully used in structural design^ and optimization 
programs^ primarily to assess the effects of the variation 
of various fundamental properties relative to the 
important physical design variables. Moreover, researchers 
have developed and applied sensitivity analysis for 
analytical model improvement and assessment of design 
trends In most cases, a predominant contributor to the 
cost and time in the optimization procedures is the 
calculation of derivatives. For this reason it is 
desirable in aerodynamic optimization to have efficient 
methods of determining the aerodynamic sensitivity 
coefficients and, wherever possible, to develop 
appropriate numerical methods for such computations 

Currently, most methods for calculating transonic 
aerodynamic sensitivity coefficients are based upon the 
finite difference approximation to the derivatives. In 
this approach, a design variable is perturbed from its 
previous value, a new complete solution is obtained, and 
the differences between the new and the old solutions are 
used to obtain the sensitivity coefficients. 

This direct, or brute force, technique has the 
disadvantage of being potentially very computer intensive, 
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especially if the governing equations are expensive to 
solve. Accordingly, the need to eliminate these costly and 
repetitive analyses is the primary motivation for the 
development of alternative efficient computational methods 
to determine the aerodynamic sensitivity coefficients 
In steady-state transonic flow problems, implicit 
approximate- factorization (AF) algorithms® have been used 
successfully and efficiently to solve the nonlinear two- 
dimensional transonic small-disturbance equation. This 
governing equation permits capturing of important physical 
phenomena while being easy to handle from a coding 
standpoint. For this reason it would seem desirable to 
initiate the study of the quas i - analytical method in the 
transonic regime using the transonic small perturbation 
equation as a practical object of this investigation 
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PROBLEM STATEMENT 

Based on the foregoing discussion, the current problem 
is formulated starting from the generic quasi-analytical 
approach and manipulated according to the rules given in 
Appendix A of Ref 1 which is reproduced in this thesis 
because of its direct significance to the derivation of 
the general sensitivity equation 

In this study, the general sensitivity equation is 
applied to the residual expression (R) of the transonic 
small perturbation equation. As mentioned earlier, this 
expression is chosen because of its simplicity as well as 
its adequate description of the nonlinear phenomena 
occurring in the transonic regime Although this 
expression is nonlinear in the perturbation potential (p), 
the general sensitivity equation, Eq (1), is linear with 
respect to the unknown sensitivity (d<p/dX D^), (see 
Appendix A) . 

It is to be noticed that the practical implementation 
of the above step is not achieved until the residual 
expression is approximated on a finite domain and the 
mathematical form of the problem rendered to that of one 
in linear algebra. This discretization process will be 
explained in a later section. 



Following the previous formulation, the . s ens i t iv i ty 


equation as applied to the residual expression of the 
transonic small perturbation equation is given by, 

( 1 ) 


where 


<3R 


' fl_. ’ 


' 3R 

dtp 


. axDi . 


. 3X Di . 


R - (B]_ + B2<p x ) Vxx + “Pyy " ® 

B — 1 — Maj 2 

B 2 - - (7+1) hao 2 

for air 7 - const. - 1.4 
<p - <p (x,y,XD) 

XD - set of design variables 
XD^ ■ i c ^ design variable 


( 2 ) 


(3) 


subject to, 

Airfoil Boundary Condition 

r d y 1 


<Py(Xb.O) - 


- F ( X , XD ) 


dx 


Infinity Boundary Condition 

subsonic : <p a , - r$/(2ir), ff 
supersonic : 9 ® - 0 , ff 

<Px “ 0 > 6 


Kutta Condition 


(4) 

njr /2 , n- 0 , 1 , 2 , 3 , 4 
njr /2 , n- 1 , 2 , 3 
njr/2 , n-0 ,4 ( 5 ) 


AP - 0 ( T - A<p - const.), xjg < x < ® 


( 6 ) 



DESIGN VARIABLES AND AIRFOIL DEFINITION 


Equation (I) is discretized into a system of linear 
equations to be solved for the unknown sensitivity vector 
The solution of this system is obtained, as explained in 
the following section, by using a Gauss-Seidel iterative 
procedure which utilizes the sparsity pattern 
characterizing the coefficient matrix (3R /d<p). An 
advantage of using this scheme is that several unknown 
vectors can be obtained simultaneously, each vector 
representing the sensitivity of the potential ( <p ) with 
respect to some design variable XD^. 

At this stage, it is convenient to define the vector 
of design variables 

XD - { XD 1( XD 2 XD n > (7) 

and to exactly determine which variables influence the 
solution of Eq.(2). In doing so, the relation between the 
sensitivity coefficients corresponding to these variables 
and the form of the optimization algorithm that utilizes 
this information needs to be considered. 

For the transonic flow problem, an appropriate choice 
of the first design variable is the free stream Mach 
Number (M*,) . This variable appears in the governing Eq.(2) 
and has an important influence on the character of the 
equation via its influence on local Mach number ( for M<1 , 



the equation is elliptic, for M>1, the equation is 
hyperbolic ) and thus on the nature of the solution For 
this reason, it is desirable to have Me as one of the 
design variables. 

Next, it is appropriate to examine the boundary 
condition given by Eq.(5). In the transonic small 
perturbation formulation, the angle of attack (a) enters 
the problem through the boundary condition and thus, 

dy 

F u “ — - y u ' - « (8) 

1 L dx Jb 1 

For simplicity, the function (F) should be easily 
differentiable with respect to the design variables 
defining the airfoil geometry. This desirable feature is 
explained later on and has to do with the computation of 
the right hand side term of the sensitivity equation. 
Therefore, it would seem plausible to have a simple 
analytical expression for modeling the upper and lower 
surfaces of the airfoil. 

For the present studies, it was decided to limit 
consideration to two basic airfoil sections, namely 
parabolic - arc sections, and the NACA four-digit sections 
These families of wing sections are obtained by combining 
a mean line and a thickness distribution^. The resultant 
expressions possess the necessary features that suit the 



problem, mainly the concise description of the airfoil 
surfaces in terms of several geometric design variables 
The expressions are as follows • 

For p ar abo 1 1 c - ar c sections 

C(2Lx-x 2 )/L 2 ± 2Tx(l-x), x < L 

C[ (l-2L)+2Lx-x 2 ] / ( 1— L) 2 ± 2Tx ( 1-x ) , x > L 

(9) 

For NACA four-digit sections 
r C(2 Lx-x 2 )/L 2 

± 5 T ( 0 29697X-0 . 126x-0 . 3516x 2 + 0 . 2843 x 3 -0 . 1015x 4 ) 

y u “ ' x < L 

1 

C[ ( 1-2L) +2Lx-x 2 ] / ( 1-L) 2 

± 5T(0 . 29697x- 0 12 6x-0 . 3 5 1 6x 2 +0 . 2 84 3x 3 -0 1015x 4 ) 

x > L 

(10) 

where 


y r 


C ■ Maximum ordinate of mean line (camber) 

L ■ Chordwise location of maximum ordinate of camber 
T ■ Maximum thickness 


Each of the quantaties C, L, and T is expressed as a 
fraction of the chord (e.g. if T is 6% chord then 
T - 0.06). Differentiating Eqs.(9) and (10) with respect 


to x and substituting the results into Eq.(8) yields . 



11 


For par abo 1 ic - ar c sections 
F u> i - 2 C ( L-x ) /LL - a 

± 2T(l-2x) (11) 

For NACA four-digit sections 
F U> 1 - 2 C ( L-x ) /LL - a 

± 5T(0 . 14845/Vx-O . 126-0 . 7032x+0 8529x 2 -0 406x 3 ) 

( 12 ) 

where 

L forward of maximum ordinate of camber, x < L 

LL “ • 

(1-L) aft of maximum ordinate of camber, x > L 

(13) 

Eqs (11) and (12) are simple analytical expressions in 
terms of the four variables T, L, C, and oc . Thus, 

XD - ( T, M a , cc, L, C } (14) 

represents the complete set of design variables that 
define the present two-dimensional airfoil sensitivity 
problem. Notice that these variables are completely 
uncoupled and, thus the sensitivity equation can be solved 


O 

independently with respect to each variable 0 . 



12 


MATHEMATICAL TREATMENT AND SOLUTION PROCEDURE 
Problem Discretization 

Equation (1) resembles the general sensitivity 
equation (see Appendix A) applied to the residual R 
instead of the function F, <p instead of y, and XD^ instead 
of x. Now, in order to solve the problem numerically, 

Eq (2) is formulated computationally on a finite domain 
This transformation is achieved by using a stretched 
Cartesian grid that maps the infinite physical domain onto 
a finite computational grid, Fig.(l). In this study, the 
grid used is based upon a hyperbolic tangent 
transformation that places the outer boundaries at 
infinity. Accordingly, the computational variables used 


are given by, 

£ - tanh A 2 x (15) 

q - tanh A^y (16) 

or , 

x - { ln[(l+0/(l-e)] ) / (2A 2 ) (17) 

y - { In [ ( l + r 7 ) / ( l-» 7 ) ] ) / (2A X ) (18) 

In terms of the grid nodes (i,j), we have, 

£ - - 1 + (i-l)Af (19) 

*?--! + (j-l)A* (20) 


Equations (15)-(20) are the equations that govern the 
hyperbolic tangent transformation. 
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Fig.l Computational and Physical Domains 
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In addition, the stretching functions are defined as, 


f - (dl/dx) - A 2 (l-£ ) 

(21) 

g - (dr?/ dy ) - Ai(1-»7 2 ) 

(22) 

that. 


cp x - fcp^ 

(23) 

<Py - g <p n 

(24) 

^xx - f (f«P^)f 

(25) 

<Pyy “ S ( g 'V 3 77 ) rj 

(26) 


Solution about a Fixed Design Point 

Substituting from Eqs . ( 23 ) - ( 26 ) into Eq.(2), yields 
the transformed residual expression, 

R - [ B 1 + B 2 f ] f(f + gigVrj),) - 0 (27) 

This equation is solved numerically by an approximate 
factorization scheme^ in which the objective is to force 
the residual to zero at each point of the computational 
domain. In finite difference form, Eq (27) can be written 
as , 

R i,j " I B 1 + B 2(‘Pi + l,j -< Pi-l, j )/(2A£) ] fi/Af 2 
["i, j f i+W(<Pi + l, j-'Pi, j ) ~ C 2 *'!, j-Dfi-h 
(<Pi, j-<Pi-l, j) “ (l-‘ / i,j) f i-3/2(‘Pi-l,j-'Pi-2,j)] 
+ (gj+H< < Pi , j+l-^i , j ) “ gj-4<*i, j“*i, j-l) J g j /Ar? 2 

(28) 
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where 

1 / 1 j - 1 if point (i,j) is subsonic 

v i j — 0 if point (i,j) is supersonic 

Eq.(28) is the discretized form of the residual at a 
general point (i,j) in terms of <p values at surrounding 
points. Consequently, R at i,j can be viewed as a function 
of the i p values at neighboring points, and, therefore, the 
differentiation of the residual expression is straight 
forward. 

Differentiation of the Residual 
Rearranging Eq.(28) yields 

R i,j “ c l c f’i , j + c 2<Pi + l , j<Pi-l , j + C3* i + 1 ,j*i,j 

+ c 4‘Pi — 1 , j ^i , j + c 5<Pi + l , jVi-2 , j + c 6<Pi-l , j^i-2 , j 
+ cy'Pi-l.j 2 + c 8<Pi + l,j 2 + c 9 <p i + li j + c lO^i-l , j 
+ cil'Pi.j+l + c 1 2 ^ i , j -1 + c 13<Pi-2,j ( 29 > 

where the coefficients C]_, c 2 , •• . c^3 are given by 

ci - ( ~ gj (gj+h + gj-4)/Ar7 2 

- fiBlt^i,jfi + h+(2^ i , J -l)f i _h]/A i e 2 ) 

c 2 - ( - n t> jB 2 fi 2 fi + h/(2A^ 3 ) 

+ fi 2 B 2 [ (2i/ lt j-l)fi_ W — (1— t / i>j )f i _ 3/2 ]/(2A^ 3 ) ) 
c 3 “ ( “ fi 2 B2[i / i, jfi+h+(2^i ) j-l)fi_h]/(2A? 3 ) ) 

c 4 - { + f i 2 B 2 [ Vi ( j f i + h + (2»'i ( j -1) f ] /( 2A? 3 ) } 

c 5 - ( + (l-^i , j )B 2 fi 2 fi_3/ 2 /(2A£ 3 ) } 



16 


c 6 

- 

{ 

- 

C 1-^i , j )B 2 fi 2 f i -3/2/(2A^ 3 ) ) 

c 7 

- 

{ 

- 

f i 2 B 2 ( ( 2 "i , j-D f i-h“< 1-^i , j ) f i-3/2 ]/( 3 ) ) 

c 8 

- 

{ 

+ 

v i , j B 2 f i 2 ?i + h/( 2 ^^ 3 ) ) 

c 9 

- 

f 

-♦* 

‘'i.jBlfifi + h/^^ 2 > 

c 10 

- 

( 

+ 

f i B lt ( 2l/ i, j-1) f i_h-(l-^i , j )fi-3/2]/^^ 2 ) 

C 11 

- 

{ 

+ 

gjgj +h/A»? 2 ) 

c 12 

- 

{ 

+ 

gjgj_li/Ar; 2 ) 

c 1 3 

- 

{ 

+ 

(l _l/ i , j ) B l f i f i-3/2/Af 2 ) 





(30) 

For 

a 

fixed computational grid, the coefficients given 


by Eq.(30) are functions only of B]_ and B 2 which in turn 
are functions of Mo, . This fact is used later when 
differentiating Eq.(29) with respect to in order to 
obtain the right hand side (3R/3M<a). 

At this stage, it is necessary to consider the 
treatment of various types of grid points and examine the 
effect on the general residual expression. As can been 
seen from the distribution of points on the computational 
domain, Fig. 2, several groups of points need special 
treatment as a result of applying various types of 
boundary conditions. Accordingly, it is necessary to 
revise the residual expression at these boundary points to 
include the boundary conditions. 
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In the following, each group of boundary points (see 
Fig 2) is denoted by a number and their corresponding 
revision is specified. It is to be noticed that applying 
the infinity boundary condition, Eq.(5), gives rise to two 
sets of updates for groups [5] thru [15]. The first 
corresponds to a subsonic free stream (i.e. if M«<1), 
whereas the second set corresponds to a supersonic free 
stream (i.e. if M <0 >1). 


Points 

Replace 

By 


[1] 

ILE<i<ITE , j-JB-1 

Vi, j+1 

vt.i + 

Ar?(y’ l- a )/gj+h 

[2] 

ILECiClTE, j-JB 

<°i,j-l 

Vi.j " 

AJ7(y'u- 0C )/gj-h 

[3] 

ITE:Si:SlM-2 , j-JB-1 

^i,j+l 

Vi, j+1 

- r 

[*] 

ITE<i<IM-2, j-JB 

Vi , j -1 


+ r 

For 

a subsonic free-stream, 



[5] 

i-IM-1 , JB<j <JM-2 

*1+1. j 

0 


[6] 

3<;i<IM-2, j-JM-1 

Vi, j+1 

- r/4 


[7] 

i-2 , 3<j^JM-2 

Vi—1 , j 

- r /2 


[8] 

3<i<IM-2, j-2 

”i. j-1 

- 3r/4 


[9] 

i-IM-1, 3<j <J B— 1 

^i+1. j 

- r 


[10] 

i-IM-1, j-JM-1 

•Pi+l.j 

0 



Pi, j+1 " r / 4 



[11] 

i-2, j-JM-1 

Vi. j+1 

- r/4 



Vi-1, j 

- r /2 

[12] 

i-2 , j -2 

Vi-1, j 

- r /2 



Vi,j-1 

- 3T/4 

[13] 

i-IM-1, j-2 

Vi, j-1 

- 3r/4 



Pi+1 , j 

- r 

[14] 

i-IM-1 , j-JB 

<Pi + l , j 

0 



«»i. j-1 

Vi, j-1 

[15] 

i-IM-1, j — J B— 1 

Vi+l.J 

- r 



Vi , j+1 

Vi, j+1 

For 

a supersonic free- 

scream , 


[5] 

i-IM-1, J B< j <JM— 2 

Vi+1, j 

Vi, j 

[6] 

3<i<IM-2 , j-JM-1 

Vi, j+1 

0 

[7] 

i-2, 3<j <JM-2 

Vi-1, j 

0 

[8] 

3<i<IM-2, j-2 

Vi, j-1 

0 

[9] 

i-IM-1, 3<j <J B-l 

Vi + 1 , j 

Vi , j 

[10] 

i-IM-1, j-JM-1 

Vi + 1, j 

Vi, j 



Vi, j+1 

0 

[11] 

i-2, j-JM-1 

Vi, j+l 

0 



Vi-1, j 

0 

[12] 

i-2, j-2 

Vi-1, j 

0 



Vi, j-1 

0 
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[13] 

i-IM-1 , 

j“2 

Vi, 2-1 

0 




Vi+1 , J 

Vi, j 

[14] 

i-IM-1 , 

j -JB 

^i+1, j 

Vi, j 




Vi, j-1 

■Pi, j-1 + r 

[15] 

i-IM-1 , 

j-JB-1 

Vi+l,} 

Vi, j 




Vi,j+l 

•Pi.j+l " r 


The above updates are used to modify the residual 
equation, Eq.(29), and yield a set of expressions, each 
being valid for a group of boundary points. The details of 
these operations are shown in Appendix B 

In setting up the complete quas i - analytical problem 
the circulation and its dependence upon trailing edge 
potentials must be carefully included. Since the 
circulation is determined by the difference in potentials 
at the trailing edge, 

r “ PuTE ~ <P 1TE ( 3i ) 

or, by interpolating the trailing edge values 

r - Ti [ 1.5 (<PITE-1,JB - ‘PlTE-l.JB-l) 

- 0.5 (<pite-1 , JB+1 “ PITE-1 , JB-2> ] 

+ T 2 [ 1.5 (vitE.JB “ <PITE,JB-1 ) 

- 0.5 (<PiTE,JB+l “ ‘PITE.JB-2 ) 1 ( 32 > 

where 

T 2 - [ £(x-0 5) - l(ITE-l) ] / (33) 

Ti - [ 1 - T 2 ] 


(34) 
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and since a branch cut extends from the trailing edge to 
downstream infinity, the trailing edge potentials appear 
in the residual expressions at points adjacent to the 
outer boundaries. Consequently, the resultant matrix 
(dR/dtp) , while banded, also contains many nonzero elements 
far from the central band. Notice that the presence of 
these elements greatly complicates the rapid and efficient 
solution of the sensitivity equation, Eq.(l), which will 
be explained later. 

Assembling (5R /d<p) and (9R/3XD^) 

The residual expressions obtained from the previous 
step are differentiated analytically with respect to the 
potential (<p). To be more specific, each equation is 
differentiated with respect to the potential at 
neighboring points and trailing edge points (the later 
enters as a result of the implicit nature of the 
circulation effects). These points are denoted by the 
counters (ii.jj) and are given by, 

( i , j - 1 ) , (i,J), ( i , j +1 ) , ( i - 2 , j ) , (i-l,j), (i + l,j), 
(ITE-1, JB-2) , (ITE-1 , JB-1) , (ITE-l.JB), ( ITE - 1 , JB + 1 ) , 
(ITE.JB-2), (ITE, JB-1), (ITE.JB), (ITE, JB+1). 

The end result is that the coefficient matrix 


(9Ri ( j/dfpn , j j ) is of size ( IM-2 ) * ( JM-2 ) x ( IM-2 ) * ( JM-2 ) 
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Again, the details and results of this step are shown m 
Appendix B. Once these relations are obtained, the actual 
coefficients are assembled by evaluating the appropriate 
analytical expresssions using a flowfield solution 
obtained from Eq.(2) for a given set of conditions (i e 
about a fixed design point). 

Similarly, the right hand side is evaluated by 
differentiating the analytical expressions for the 
residual (see Appendix B) with respect to each design 
var i able. 

Solution by Gauss-Seidel 

For a general (IM*JM) grid, the system given by Eq.(l) 
is of size ( IM-2 ) * ( JM-2 ) x ( IM-2 ) * ( JM-2 ) . This system is 
large, of block structure, and sparse, and, as mentioned 
earlier while banded, also contains many nonzero elements 
far from the central band As a result of this size and 
structure, and since the primary objective of this study 
is to establish the feasibility of the quas i - analy t ical 
method, it was obvious that a reasonably efficient scheme 
for solving Eq.(l) was needed and that an elimination 
technique, while straightforward, would be too time 
consuming. Consequently, the results presented in this 
thesis have been obtained using a Gauss-Seidel iterative 
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scheme^. This scheme has not been optimized for speed 
(through the choice of optimum acceleration parameters) 
but uses sparse matrix technology in processing only the 
nonzero elements; and, therefore, it is considerably 
faster and more efficient (with regard to storage 
requirements) than elimination methods. 

In handling the sparsity pattern, the symbolic 
assembly of the coefficient matrix is performed only once 
for a given grid size and given free-stream (subsonic 
versus supersonic). The resultant structure is then stored 
on a diskfile. Before the numerical part is executed, the 
symbolic information is read into the code and used 
directly to assemble the new matrix. This procedure is 
followed in order to reduce the time consumed in 
assembling the coefficient matrix. 

Once the s ens i t i vi t ies of the potentials, and thus the 
Cp distribution, to the design variables are known, the 
sensitivity of the lift coefficients to the design 
variables can be easily computed. To minimize errors, 
these coefficients are computed using 

C L - 2 T - 2 (<p u TE _< P1TE) (35) 

and henc e , 

SCl/3XD - 2 (dcp u j£/dXD i -3<p^ TE / axD i) (36) 



TEST CASES 


In this study, the quas i - analytical method has been 
used to determine the aerodynamic sensitivity coefficients 
at three freestream Mach numbers ( 2, .8, and 1 2) for two 

arbitrarily selected airfoils, each at one degree angle of 
attack. The first is a cambered parabolic arc section 
having 1% camber at 40% chord, a maximum thickness of 6% 
at 50% chord, and which is designated P1406; and the 
second is a NACA 1406 airfoil. Lift coefficients computed 
for these cases are shown in Table 1. 

In the following, two types of results will be 
presented. The first will be plots of Cp versus chord for 
the three chosen Mach numbers and two airfoil sections. 

The second will be the corresponding plots of (5Cp/3T), 
(aCp/aMa,), (acp/aoc) , (acp/ac), and (aep/aL) obtained by 
the quas i - analytical method. In addition, all of the 
figures will also contain results obtained using the 
direct (finite difference) approach in which each design 
variable was individually perturbed by a small amount, 
typically 0.001, and a new flowfield solution obtained. 
Then the sensitivities were computed using ACp/AXD . 

Finally, tables containing lift coefficients, lift 
coefficient sensitivity coefficients, and time comparisons 
are presented for all cases. 












In all subsonic cases (M® 


0 2 , 0 8 ) , an 81*20 


stretched Cartesian grid was utilized. For the supersonic 
case (Mo, - 1.2), a 41*20 grid was used. In addition, for 
these studies, the flowfield was normally computed using 
double precision arithmetic and the maximum residual 
reduced eight orders of magnitude. It was felt that this 
level of convergence was necessary in order to accurately 
evaluate sensitivity coefficients using a finite 
difference approach, although such convergence may not be 
required in the flowfield slover for the quas i - ana ly t i c a 1 
method. 

Notice also that in all cases the error tolerances 
used in the Gauss-Seidel solver for the coefficients 
involving maximum thickness, free stream Mach number, and 
location of maximum camber were l.E-06 while those on 
angle of attack and maximum camber were 1 . E-04 . 
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RESULTS AND DISCUSSION 
Accuracy of the Quas i -Analytical Method 

In order to verify the accuracy of the quasi- 
analytical method, design sensitivity coefficients for Cp 
and Cl were obtained using the finite difference 
approach. In this procedure, a single design variable was 
perturbed by typically 0.001, while all others remained 
constant, and a new flowfield solution obtained using the 
approximate factorization solver. Then values of the 
various sensitivity coefficients were obtained by finite 
differences. The results obtained in this manner for the 
pressure distributions have been shown on all figures by 
dashed lines, and in many cases the dashed lines are 
coincident with the quas i - analytical results (solid 
lines). In addition, Table 2 compares results obtained by 
the two methods at Mach numbers 0.2, 0.8, and 1.2. In most 
cases the agreement is within significantly less than 1% 

Subsonic Cases (M«, — 0.2) 

P1406 Airfoil: Initial studies concentrated on subsonic 
cases since such cases should run quickly and since at 
least approximate results would be known from thin airfoil 
theory. Figure 3 shows the pressure distribution for the 
P1406 airfoil while Figs. 4a and 4b show the sensitivity of 
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Table 2 

Accuracy of Quasi-Analy tical Method for Computing 
Lift Coefficient Sensitivity Coefficients 


P1406 Airfoil, Grid 81*20 


XDi 

M. = 

0.2 

M. = 

0.8 

M. = 

1.2 * 

QA 

FD 

QA 

FD 

QA 

FD 

T 

0.0050 

0.0053 

1 . 1208 

1 .1177 

-0.2505 

-0.2476 

M. 

0.0472 

0.0476 

1.4556 

1 . 4690 

-0.1259 

-0.1252 

a 

6.1075 

6.1383 

1 0 . 6667 

10.7742 

5.0899 

5.0920 

c 

9.9140 

9.9434 

20.5295 

20.4958 

1.3152 

1.3279 

L 

0.0692 

0.0697 

0.1647 

0.1663 

-0.1121 

-0.1114 

NACA 

1406 Airfoil, Grid 81*20 




XDi 

M- = 

0.2 

M« = 

0.8 

M_ = 

1.2 * 

QA 

FD 

QA 

FD 

QA 

FD 

T 

0.0044 

0.0044 

0.5433 

0.5231 

-0.3322 

-0.3305 

M. 

0.0466 

0.0471 

0.9891 

0.9708 

-0.0803 

-0.0802 

a 

6.1077 

6.1386 

10.3861 

10.5227 

4.5982 

4.5983 

C 

9.9086 

9.9380 

18.4786 

19.5767 

1.3495 

1.3769 

L 

0.0690 

0.0696 

0.1482 

0.1499 

-0.0304 

-0.0301 


* Executed on Grid 41*20 
QA Quasi-Analytical 


FD Finite Differece 
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the pressure to thickness for the same airfoil. As 
expected from thin airfoil theory, the upper and lower 
surface values are essentially identical and the 
difference is very small everywhere. Also shown on the 
same figure (and on subsequent figures) by the dashed line 
is the result obtained by using the finite difference 
approach; and as can be seen, the agreement betweeen the 
two approaches is excellent. 

The sensitivity of pressure to freestream Mach number 
is plotted on Figs . 5a and 5b. It is noticed that while the 
profiles for the upper and lower surfaces are similar, 
they are not equal in magnitude, indicating a nonlinear 
variation with Mach number as predicted by simple Prandtl- 
Glauret Theory. However, as indicated by the results on 
Fig. 5b, the magnitudes for this subsonic Mach number are 
very low. 

The sensitivity of the pressure coefficients to angle 
of attack are shown for this case on Figs. 6a and 6b. As 
expected from linear thin airfoil theory, the upper and 
lower surface curves are essentially equal in magnitude 
but of opposite sign. Not surprisingly, the sensitivity of 
the delta Cp variation, Fig. 6b, has the shape of the 
pressure difference curve for a flat plate at angle of 




(b) 

Fig. 5 Sensitivity of Pressure to Mach Number, 
P1406 Airfoil, - 0.2, « - 1 • 





Fig. 6 Sensitivity of Pressure to Angle of Attack, 


P1406 Airfoil, M„ - 0.2, « - 1 * 
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attack; and its magnitude, particularly near the leading 
edge is quite large 

On Figs 7a and 7b is plotted the sensitivity of the 
pressure coefficient to the amount of maximum camber 
Since camber contributes to lift, it is expected from thin 
airfoil theory that these values should be "equal but 
opposite in sign" for the upper and lower surfaces In 
addition, the pressure difference curve has the correct 
shape for that associated with a 14 mean line with the 
peak occuring at 30% chord^ and has magnitude comparable 
to those for the (3Cp/3<x) curves. 

Finally, the sensitivity of pressure to the location 
of the maximum camber point is portrayed on Figs. 8a and 
8b, and to say the least the results are interesting. 

Since maximum camber location affects the camber profile 
and hence lift, the equal and opposite behavior of the 
upper and lower surface coefficients is expected. In 
addition, the pressure difference sensitivity is primarily 
negative forward of the point of maximum camber and 
positive aft of it. This result indicates that if the 
location of maximum camber were moved rearward slightly 
(i.e. a positive AL) that lift would be decreased on the 
forward portion of the airfoil and increased on the aft 
portion of the airfoil, which is in agreement with the 
results presented in Ref. 7. 












NACA 1406 Airfoil: Figure 9 shows the pressure 
distribution for the NACA1406 airfoil at M® - 0 2, and it 
should be noted that it is different than that obtained 
for the P1406 airfoil. This difference is basically due to 
the different thickness distributions describing each 
profile. Figs. 10 to 14 show for the NACA 1406 airfoil the 
sensitivity of ACp with position for each of the five 
design variables. Since the NACA 1406 and the parabolic 
P1406 both have the same camber line and since for this 
low Mach number and thin airfoils the solutions should 
essentially be linear, the sensitivity to maximum camber 
and location of maximum camber should be essentially 
identical for the two airfoils. As can be seen by 
comparing Figs. 7 with 13 and 8 with 14, the present quasi- 
analytiacl method does indeed yield this result. Likewise, 
the sensitivity to angle of attack, Figs. 6 and 12, are 
also identical for the two airfoils. 

However, the pressure sensitivity to thickness, Fig. 

10, and freestream Mach number, Fig. 11, while very small 
in magnitude compared to the other coefficients, have a 
different chordwise variation than that for the P1406 
airfoil. The first, of course, is expected since the two 
airfoils have different thickness distr ibutions ; and the 
second is due to the fact that the two airfoils have 
























-2i 

Fig. 13 Sensitivity of Pressure to Maximum Camber, 
NACA 1406 Airfoil, M« - 0.2, <x - 1 * 







(a) 



(b) 

Fig 14 Sensitivity of Pressure to Location 
of Maximum Camber, 

NACA 1406 Airfoil, M® - 0.2, « - 1 ' 
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entirely different pressure distributions and, thus, 
sensitivity to Mach number. 

The sensitivity of the lift coefficients to the design 
variables were shown in Table 2. It should be noticed that 
for the subsonic case that the lift sensitivities for the 
two airfoils are essentially identical. Since these 
airfoils are thin and since they have the same camber 
line, such agreement should exist. 

Transonic Cases (M„ - 0.8) 

P1406 Airfoil: For this case, the cambered parabolic 
airfoil is slightly supercritical with a weak shock on the 
upper surface at about 55% chord, Fig. 15; and the lower 
surface is entirely subcritical with the minimum pressure 
point ocurring at 60% chord. As a consequence, the 
variation with chord of the sensitivity coefficients is 
considerably different than in the subsonic case. 

Figs. 16a and 16b show the sensitivity of pressure to 
the maximum thickness; and while the lower surface profile 
is similar to that obtained at subsonic conditions, the 
upper surface curve and the pressure difference 
coefficient plot show the effect of the upper surface 
shock wave. The large peak on the curves corresponds to 
the location of the shock wave and indicates that the 












shock wave location is very sensitive to maximum 
thickness Notice on Figs 16a and 16b the excellent 
agreement of the quas 1 - ana ly 1 1 c a 1 results indicated by the 
solid lines with those obtained using the f ini t e - d i f f e r e c e 
approach (dashed lines) 

The results for (3Cp/3M „), which are shown on Figs 17a 
and 17b, are similar The lower surface curve is typical 
of a subsonic flow, while the upper surface and the 
pressure difference coefficients reflect the presence 
of the upper surface shock wave Similar comments can be 
made for the remaining design variable coeficients, which 
are plotted on Figs 18, 19, and 20. 

Examination of the curves in the vicinity of the shock 
wave location indicates that the pressure sensitivity and 
indirectly the shock wave location is about equally 
influenced by the maximum thickness, freestream Mach 
number, and angle of attack. 

However, in comparison it is relatively insensitive to 
location of maximum camber; but, perhaps surprisingly so, 
the pressure is twice as sensitive to the amount of 
maximum camber as it is to the other design variables. 

NACA 1406: At M*, - 0 8, the flow about the NACA 1406 
airfoil has a strong shock at 40% chord, Fig. 21. As a 
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(b) 


Fig. 18 Sensitivity of Pressure to Angle of Attack, 


P1406 Airfoil, M* 


0.8, oc 


1 
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result, the pressure design sensitivity coefficients, 
which are shown on Figs. 22-26, have large peaks at the 
shock location. In addition, while the curves are similar 
to those obtained for the P1406 airfoil, they differ in 
details and some cases in magnitude. In particular, the 
peak value at the shock in the 3Cp/3<x curve for the NACA 
1406 is significantly higher than that for the P1406 
airfoil Further, for the subsonic "linear" flow situation 
the sensitivity coefficients for angle of attack, maximum 
camber and location of maximum camber were identical for 
the two airfoils. However, at transonic conditions, the 
flow is highly nonlinear and the corresponding curves for 
the two airfoils are significantly different. Again the 
reasonable agreement between the quas i - analytical results 
(solid lines) and the finite difference results (dashed 
lines) is evident on the figures. 

The sensitivity of lift to the design variables was 
also shown in Table 2 for both airfoils. While the values 
for the various design variables are similar in magnitude 
for the two airfoils, there are some significant 
differences. For example, the coefficients for maximum 
thickness differ by a factor of two between the two 
airfoils and those for Mach number differ by about fifty 
percent. Also, it should be noticed that for both airfoils 










(b) 


Fig. 23 Sensitivity of Pressure to Mach Number, 
NACA 1406 Airfoil, M* - 0.8, ce - 1 ' 





(a) 



(b) 


Fig. 24 Sensitivity of Pressure to Angle of Attack, 
NACA 1406 Airfoil, M« - 0.8, « - 1 * 
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Fig. 25 Sensitivity of Pressure to Maximum Camber, 
NACA 1406 Airfoil, - 0.8, oc - 1 * 







Fig. 26 Sensitivity of Pressure to Location 
of Maximum Camber, 

NACA 1406 Airfoil, M„ - 0.8, « - 1 * 





the lift is most sensitive to angle of attack and to 
maximum camber. In addition, Fig. 25 shows a discrepancy 
between the results obtained by the direct approach and 
those obtained thru the quas i - analytical method This 
discrepancy could be due to either the magnitude of the 
perturbation used to compute the sensitivities by the 
direct approach, or the tolerance values used in the 
stopping criteria of the Gauss-Seidel procedure. 

C. Supersonic Cases (M« - 1.2) 

In order to investigate the applicability of the 
quas i - analytical method at supersonic freestream Mach 
numbers, solutions were obtained for the P1406 and the 
NACA 1406 airfoils at Mach 1.2. At this condition, the 
flow is transonic in that the bow shock is detached, and 
there is a region of subsonic flow extending to 
approximately the quarter chord, Figs. 27 and 33. Figures 
28-32 and 34-38 show the pressure sensitivities for these 
cases, and Table 2 listed the lift sensitivities. 

As mentioned earlier, notice that for these cases the 
solutions presented are for the 41*20 medium grid. 

Examination of the plots shows that the pressure 
sensitivity coefficients have different trends and 
magnitudes from those computed for subsonic freestream 









































a Moo 
- 1.0 



(a) 





- 2 ) 

Fig 36 Sensitivity of Pressure to Angle of Attack 
NACA 1406 Airfoil, M* - 1.2, ce - 1 * 
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Fig. 37 Sensitivity of Pressure to Maximum Camber, 
NACA 1406 Airfoil, M* - 1.2, « - 1 0 














72 


supercritical conditions, Figs. 16-20 and 22-26, and that 
they are approaching the form expected from supersonic 
linear theory. These changes are particularly evident in 
the lift derivatives presented in Table 2. Notice that the 
derivatives with respect to the design variables maximum 
thickness, Mach number, and location of maximum camber 
have switched sign. In addition, as expected from linear 
theory, the influence of camber on lift has decreased 
significantly; and at M«, - 1.2 is only about 25% of the 
angle of attack effect as compared to a factor of about 
two at H b - 0.8. 

Time Comparisons 

Obviously, in the development of the quas i - analytical 
method it was hoped that not only would this approach 
yield accurate values for the aerodynamic sensitivity 
coefficients but also that it would be more efficient than 
the brute force finite difference approach. Tables 3 and 4 
present some comparisons concerning the amount of 
computational effort required to obtain solutions by the 
two approaches . 

In comparing the values, several items should be kept 
in mind. First, it has been assumed that the finite 
difference approach will require six independent 



Time Comparisons for Obtai 


for Five Design Va 
P1406 airfoil 


Method 


Flowfield Solver 

Finite Difference Approch 
(6 Flowfield Solutions) 

Quasi-Ana lytical Approach 
(1 Flowfield Solution plus 
Sensitivity Coefficient 
Solution via Gauss-Seidel) 

Ratio QA / FD 


NACA1406 airfoil 


Method 


Flowfield Solver 

Finite Difference Approch 
(6 Flowfield Solutions) 

Quasi-Analytical Approach 
(1 Flowfield Solution plus 
Sensitivity Coefficient 
Solution via Gauss-Seidel ) 

Ratio QA / FD 


Note: All times are normalize 


at Mach 0.2. 




















Table 4 

Time Comparisons for Obtaining Sensitivity Coefficients 
for Five Design Variables, Grid 81*20 
P1406 airfoil 


Method 

M. = 0.2 

M. = 0.8 

Flowfield Solver 

1.00 

1.93 

Finite Difference Approch 
(6 Flowfield Solutions) 

6.00 

11.58 

Quasi-Analy tical Approach 
(1 Flowfield Solution plus 
Sensitivity Coefficient 
Solution via Sauss-Seide 1 ) 

5.56 

14 .65 

Ratio QA / FD 

0.93 

1.27 


NACA1406 airfoil 


Method 

3 

1 

II 

O 

K) 

M. = 0.8 

Flowfield Solver 

1.00 

2.46 

Finite Difference Approch 
(6 Flowfield Solutions) 

6.00 

14.76 

Quasi-Analytical Approach 
(1 Flowfield Solution plus 
Sensitivity Coefficient 
Solution via Gauss-Seidel ) 

7.39 

26.29 

Ratio QA / FD 

1.23 

1.78 


e: All times are normalized by the flowfield 
solver at Mach 0.2. 
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solutions. In practice it might be possible to start each 
finite difference solution from a previous solution and, 
thus, decrease the time to convergence. However, to be 
accurate, the finite difference approach will probably 
require double precision and will have to be extremely 
well converged (i.e.l.E-08) Nevertheless, the values for 
the finite difference approach probably should be viewed 
as maximum values. 

Second, the Gauss-Seidel method for obtaining the 
sensitivity coefficients has not been optimized and may 
not even be an efficient method; and the flowfield 
solution required for the quas i - analytical approach may 
not need double precision and may not have to be as 
tightly converged. Thus the values shown for the quasi- 
analytical approach should also be viewed as maximum 
values . 

In spite of these limitations, the results do 
indicate, for the 41*20 grid, that the quas i - analytical 
method is at transonic conditions potentially more 
computationally efficient than the brute force finite 
difference approach. 

As for the fine grid (81*20), it is seen that the 
finite differece method is generally more efficient than 
the quas i - analytical method. This result implies the need 



to improve the efficiency of the sensitivity coefficient 
solver, which could be achieved thru the use of an optimum 
acceleration procedure Another alternative explanation 
could be that the Gauss-Seidel solver might not be an 
efficient algorithm in the case of fine grids. This 
possibility will be discussed in the following section. 



CONCLUSIONS AND RECOMMENDATIONS 


Based upon Che investigations and results described i 

the previous sections, it is concluded that 

(1) The Quas i -Analytical Method is feasible as applied to 
the transonic small perturbation residual expression 

(2) The results obtained from the quas i - analytical method 
are almost identical to those obtained by the brute 
force technique (finite difference). 

(3) The Gauss-Seidel iterative procedure used in solving 
the sensitivity equation is potentially more 
computationally efficient than the brute force direct 
approach for the 41*20 medium grid. 

(4) For fine grids (81*20), the Gauss-Seidel solver is no 
as competi t i ve . 

It is recommended to 

(1) Conduct more studies using refined grids at a variety 
of free stream Mach numbers. 

(2) Consider new solution schemes for the quasi- 
analytical equation, Eq.(l) 

(3) Examine the concept of only including part of the 
flowfield in attempt to reduce the size of system (1) 

(4) Extend the calculation of the lift coefficient 
sensitivity derivates to that of the moment 
coefficient sensitivity derivates. 
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It is also recommended that further studies be carried 
out to examine the structure and properties of the 
coefficient matrix for special cases, as for example, that 
of the non-lifting problem (°c - 0, C - 0), and accordingly 
determine the feasibilty and efficiency of new fast 
solvers that exploit this structure. In addition, other 
iterative tridiagonal solvers need to be investigated 
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APPENDIX A 

GENERAL EQUATION FOR SENSITIVITY 
This Appendix^- is a self-contained tutorial on 
sensitivity analysis arising in a generic problem whose 
governing equations are given. Let 

F(y , x) - 0 (Al) 

represent governing equations of a problem in which y is 
unknown to be obtained by solving Eq.(Al), and x are given 
constants. The quantities y and x may be vectors, and F 
may be a vector of functions If y is a vector, Eq (Al) 
implies a set of equations whose number is equal to the 
length of vector y; however, the x vector may be shorter 
than y. Existence of the solution of Eq (Al) makes, 
implicitly, y - f(x). The functions F may be anything 
computable : linear algebraical equations, partial 

differential equations, integral equations, or integral- 
differential equations, transcendental functions, etc It 
may be nonlinear, and may require an iterative method for 
solution of Eq.(Al). 

If Eq.(Al) governs a physical system being designed, 
then the designer wants to know not only the y for a given 
x, but also the sensitivity of y to those x-quantities 
that he controls as design variables. For instance, F(y,x) 
might be the Euler equations from which to compute y - the 
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pressure distribution on a body in airflow, and x might be 
the body geometry variables The designer of the body 
shape needs to know 3y/3x 

One way to obtain 3y/3x is by finite differences This 
requires solving Eq (Al) for a given x to obtain y Then 
assume, for one element of x, a perturbation, x - x + Ax, 
and repeat solution of Eq (Al) to get y + Ay. Then, an 
approximation to 3y/3x is 

3y/3x - Ay/Ax (A2) 

This operation must be repeated for all x-quantities of 
interest and may be prohibitively computer - intens ive , if 
Eq.(Al) is expensive to solve. In addition, the accuracy 
of 3y/3x will depend on the proper choice of Ax. 

An alternative is the quas i - analytical approach. It is 
called "quasi-" because the y(x) is known only 
numerically. However, for Ax, it must be true that 

F(y+Ay, x+Ax) - 0 (A3) 

In other words, increase of x must be compensated for by a 
change in y to preserve the zero value of F. Hence, 
recognizing that the total derivative of f with respect to 
x is according to the textbook rules of differentiation 
for implicit functions 

dF/dx - 3 F/3x + (3F/3y) (3y/3x) (A4) 
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Eq (A2) will be satisfied if 

(dF/dx) Ax - 0 (A5) 

Substituting Eq.(A4) into (A5), and rearranging, yields 

( 3 F/3y ) ( 3y/3x) (3F/3x) (A6) 

Eq . ( A6 ) is a general sensitivity equation in which the 
desired sensitivity appears directly as the unknown 
(3y/3x). For a vector y of length n, the term (3F/3y) is a 
matrix n*n whose each column is a vector of gradients with 
respect to y (a Jacobian matrix), the term (3y/3x) is a 
vector of unknown derivatives of y with respect to one 
particular x variable, and the term (3F/3x) is a vector of 
derivatives with respect to the same particular variable 
x. Computation of the derivatives of y with respect to 
several variables x requires solutions of Eq.(A6) with 
many right hand sides, one per each variable x. Since the 
Jacobian matrix remains the same for all variables x, a 
solution algorithm arranged so as to factor the matrix 
only once will be preferred for computational economy. 

It is important that Eq.(A6) is simply a set of 
linear, algebraical equations even though Eq.(Al) may be 
far more complicated than that. The terms (3F/3y) and 
(3F/3x) may still not be obtainable analytically. If so, 
they can be computed by finite difference, i.e , assuming 
perturbation x - x + Ax and y - y + Ay for each element of 
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x and each element of y separately, and substituting into 
Eq (Al) , one obtains the respective AF values (upon 
substitution of x + Ax, or y + Ay, F in Eq (Al) is no 
longer equal to zero, it becomes AF) from which the terms 
(3F/3y) and (3F/3x) can be computed as in Eq (A2) 

Computation of the terms (3F/3y) and (3F/3x) by finite 
difference is accomplished by repetitive evaluations of 
F(y,x) for known y and x, as opposed to repetitive 
solutions of Eq (Al) for unknown y required by Eq (A2) 
Hence, the quas 1 - ana ly t i c a 1 approach is inherently less 
computer intensitive than the finite difference procedure 
based on Eq . ( A2 ) . 



APPENDIX B 


RESIDUAL EXPRESSIONS AND DERIVATIVES 
The residual expression at a general point (i,j) is 
given by , 

R i,j “ c l*i,J + c 2^i + l , j Pi-1 , j + CB^i + l , j<Pi , j 

+ c 4'Pi-l 1 jVi i j + C5<p i + i ( j<Pi_ 2l j + c 6^i-l , j^i-2 , j 
+ c 7<Pi-l,j 2 + c 8*Pi + l , j * + cg'Pi + l.j + c 10'Pi-l , j 
+ c 1 l^i , j +1 + c 12‘Pi,j-l + c 13*i-2,j 
For points designated [l]-[4], we have, 

[1] R lt j - (ci + cix)«pi f j + c 2 < i p i + l , j •Pi-l , j + c 3<Pi + l , j*i , j 

+ c 4<Pi-l , j Vi , j + c 5Vi + 1 ( j Vi _ 2 , j + C6<Pi-l, j Vi — 2 . j 

+ c 7 < Pi-l , j 2 + cSVi + l.j 2 + c 9 «p i + 1( j + c 10 dPi_i ( j 
+ c 12<Pi,j-l + c 13*>i-2,j + cnAr ? F 1 /g J+ 4 

[2] R i( j - (ci + c 12 )?>i , j + c 2 9i + l , j<Pi-l , j + c 3 ¥> i + i , j<Pi , j 

+ c 4^i-l, j*i, j + c 5^i + l , j ^i-2 , j + c 6<Pi-l , j<Pi-2 , j 

+ c 7 < Pi-l,j 2 + c S'Pi + l , j 2 + c 9<Pi + l,j + c 10«Pi-l,j 
+ c ll < Pi , j +1 + c 13Vi-2,j " c 12 A, 7 F u/gj-h 

[3] R it j - c 1( p i( j + c 2 «p i + 1> j<Pi-l, j + c 3 <p i + 1 f j<Pi f j 

+ c 4?i-l , j?i , j + c 5?i+l , j ^i-2 , j + c 6<Pi-l , j'Pi-2 , j 
+ c 7?i-l,j 2 + c 8‘i p i + l,j 2 + c 9<Pi + l,j + c 10<Pi-l,j 
+ c ll^i , j +1 + c 12?i,j-l + c 13?i-2,j " c ll r 

[4] R i( j - c 1<Pi> j + c 2 ^ i + 1 _ jiPi-x t j + c 3 ^ i + 1( jVi,j 

+ c 4<Pi-l , j<Pi , j + c 5‘Fi + l , j ‘Pi^ ,j + c 6‘Pi-l,j < Pi-2,j 
+ cy^x.xj 2 + cg^i + XJ 2 + c 9<Fi + l,j + c 10<Pi-l,j 
c ll ( Pi , j +1 + c 12?i,j-l + c 13<Pi-2,j + c 12 r 


+ 
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For a subsonic free stream, 

[5] R i,j " c l'Pi , j + c 4 < Pi-l , j^i , j + c 6?i-l , j^i-2 , j 

+ c 7<Pi-l,j 2 + c 10 < Pi-l , j + c 1 l^i , j +1 
+ c 12^i,j-l + c 13<Pi-2,j 

[6] R if j - c 1( p 1( j + c 2 <Pi + i t j<Pi-i,j + c 3^i + l , j<Pi , j 

+ c 4<Pi-i , jVi , j + c 5 p i + 1 , jiPi-2 , j + c 6 | Pi-l , j “Pi-2 , j 
+ C 7<P i-1 , j 2 + c 8^i + l,j 2 + c 9<Pi + l,j + c 10^i-l,j 
+ c u (-r/4) + c 12 «Pi , j _i + c 13t p i _ 2> j 
( 7 ) R i,j “ c lVi,j + c 2 «Pi + l , j (-T/2) + c 3 <p i + 1 f j^i , j 

+ c 4 (-r/2)vi, j + c 5 ( Pi + l , j ^i -2 , j + c 6 (-r/ 2 )<p i _ 2 _ j 
+ c 7 (-r/ 2 ) 2 + cs^i + ij 2 + c 9 ( Fi + l,j + ci 0 (- r /2) 

+ c ll<Pi,j+l + c 12*>i,j-l + c 13Vi _ 2> j 
[8] R i( j - c 1<Pi| j + c 2 <p i + 1 f j<Pi_i , j + c 3 «p i + i f j<pi ( j 

+ c 4<?i-l , jVi , j + c 5 ^ i + i > j^i_2 , j + c 6 ( Pi-l , j Vi-2 , j 
2 2 

+ C7<Pi_ lfJ + c 8 <p i + 1> j + C9<p i + 1 j + c 10 «pi_i ( j 

+ c 1 l^i , j +1 + c 12 (-3r/4) + ci2Vi-2,2 
[ 9 1 R i , j - c l < Pi , j + c 2 (- r )<Pi-l,j + csC-Dvi.j 

+ C4<Fi-l , j<Pi , j + c 5 (-r)«Pi_ 2 , j + c 6 <Pi_i ( j^i-2 , j 
2 2 

+ c 7 <Pi-l,j + c 8 (-r) + cg(-r) + ciovi-l.j 

+ c 1 l^i , j +1 + c 12<Pi,j-l + c 13 < Pi-2 , j 
[10] R i( j - c 1( p i( j + c 4 p i _ li j9 i(J + c 6V >i_i , j<?i-2 , j 
+ c 7ipi _ lt j 2 + c 10 «iPi-l,j + c 11 (-r/4) 


c 12?i,j-l + c 13?i-2,j 
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[11] R i , 


[12] R i, 


[13] Ri, 


[1^] Ri, 


[15] Ri, 


j “ c l^i,j + c 2^i + l , j ( -r /2 ) + cS'Pi + l , j ^i , J 
+ c 4 (-r/2) ¥ » itJ + c 5 <p i + i t 2 <Pi -2 , j 
+ eg ( — r / 2 ) i _ 2 , j + c 7 (“T/2 ) + c 8 l Pi + l , j 
+ C 9 <pi + i ( j + c 10 (-r/2) + cn(-r/4) 

+ c 12 <Pi,j-l + c 13‘Pi-2 , j 
j “ ci^ij + c 2‘Pi + l , j (“ r / 2 ) + c 3'Pi + l, j'Pi.j 
+ C4(-r/2)*i, j + c 5 <p i + i ( jVi_2 , j 
+ eg (-T/2 )<Pi-2 , j + c 7 (“ r /2 ) + C8«Pi + l,j 

+ c 9 <p i + lj j + c 10 (-r/2) + cuvi.j+l 
+ c 12 (-3r/4) + c 13 <Pi_ 2 ,j 
j “ c l^i,j + c 2 (-r)vi-l , j + c 3 (“R) < Pi , j 

+ c 4 .pi_i t j <Pi 1 j + C 5 < — r ) <p i_ 2 , j + C6“Pi-l , j'Pi- 2 , J 

+ c 7<Pi-i,j 2 + cg(-r) + C 9 (-r) + cxo'Pi-l.j 
+ c ll<Pi,j + l + c 12 (-3r/4) + c 13 <Pi_ 2> J 
j “ c l^i , j + c 4?i-l, j<Pi, j + c 6^i-l , j'Pi-2 , j 
+ c 7*Pi-l , j 2 + c 10^i-l , j + c ll l Pi , j + 1 
+ c 12 <p i( j_i + c 13 (p i _ 2t j + c 12 r 

j “ c l < Pi , j + c 2 ( -r ) ‘P i- 1 , j + c 3(~ r )'Pi,j 

+ c 4 (Pi_x f j <Pi t j + c 5 (-r)«pi -2 , j + c 6*Pi-l , j Pi-2 , j 
+ C7<pi_i f j 2 + cg(-T) + C 9 (-r) + cio'Pi-l.j 
+ c ll^i , j + 1 + c 12<PiJ-l + c 13'Pi-2,j “ c ll r 



For a supersonic free stream, 


2 


[5] Ri, 


[6] R it 


[7] R it 


[8] Ri, 


[9] Ri, 


[10] Ri 


j ” c l^i , j + c 2<Pi , j ^i-1 , j + c 3^i,j 

+ c 4 <p i _ 1 , jVi , j + c 5<Pi, jVi-2, j + c 6<Pi-l , j^i-2 , j 

+ c 7<Pi-l,j 2 + c 8<Pi,j 2 + c 9<Pi,j + c 10<Pi-l,j 

+ c 1 l^i , j +1 + c 12 < Pi , j -1 + c 13<Pi-2,j 

j ” c l^i , j + c 2 < Pi + l , j*i-l , j + c 3 ( ( , i + l , , 3 

+ c 4<Pi-i, j<Pi, j + c 5 <p i + 1 , j<Pi_ 2 , j + C6<Pi- lp j*i-2,j 
2 2 

+ c 7^i — 1 , j + C8<Pi + l,j + c 9 <p i + 1 ,j + c 10 <Pi-l,j 

+ c 12 <Pi,j_i + c 13 <pi_ 2 ,j 

j - c l'Pi , j + c 3^i + l , j'Pi , j + c 5<?i + l , j^i-2 , j 

+ c 8 l Pi + l , j 2 + c 9<Pi + l,j + cil^ij+l 

+ c 12^i , j -1 + c 13<Pi-2,j 

j “ c l^i , j + c 2<Pi + l , j Vi — l , j + c 3 <p i + 1 , jVi , j 

+ c 4<Pi-l , j<Pi , j + c 5 <p i + 1 , j<Pi_ 2 , j + c 6 <pi_i, jPi_ 2 ,j 
2 2 

+ cy^i.ij + c 8 <p i + 1>J + cgvi + l.j + c 10 <Pi-l,j 

+ c ll < Pi , J +1 + c 13 ( Pi-2,j 

j “ c 1<P i , J + c 2^i, j^i-l.j + £3^,^ 

+ c 4¥ > i _ 1 , j^i , j + c 5Pi , j SPi — 2 , j + c 6«>i-l, j?i-2, j 
+ c 7 <pi_ 1 ( j 2 + c 8 <Pi,j 2 + c 9 <p ij + ciovi-i.j 
+ c ll'Pi , j + 1 + c 12<Pi,j-l + c 13<Pi-2,j 
,j “ ciVi.j + c 2 < Pi , j *Pi-l , j + c 3<*>i,j 2 

+ c 4?i — 1, jv»i, j + c 5?i, j?i-2 , j + c 6^i-l, j«Pi-2 , j 
+ cy^i.ij 2 + c 8 <pi,j 2 + c 9 «p i( j + c 10 <Pi_i,j 


+ c 12 ^i,j_l + c 13 <pi_ 2 ,j 
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[ 1 . 1 1 R i,j “ ci'Pi.j + c 3<Pi + l , j<Pi , j + c 5<Pi + l , j^i-2 , j 

+ c 8^i + l , j 2 + c 9<Pi + l,j + c 12 ( Pi 1 j-l + c 1 3^1-2 , j 
t 12 ] R i , j ” c l«Pi,j + c 3<Pi + l, J*i, j + c 5 <p i + 1 , j<Pi-2 , j 

+ c 8<Pi + l,j 2 + c 9^i + l,j + c ll‘Pi,j+l + c 13‘Pi-2 , j 
[ 13 1 R i,j ” c l<Pi ( j + c 2 ( Pi , j SP i — 1 , j + c 3*>i,j 2 

+ c 4Vi-l , j<Pi , j + =5^1 , j-Pi-2 , j + c 6<Pi-l , j^i-2 , j 
+ C 7<P i - 1 , j 2 + c 8 sp i , j + c 9'Pi,j + c 10 < Pi-l , j 
+ c ll ( Pi , j +1 + c 13<Pi-2,j 

[14] R it j - c 1<Pi| j + c 2 «Pi , j<Pi-l , j + c 3 <p i( j 2 

+ c U<PL-l , jVi , j + 0 5 'P i. , j Vt-2 , j + c 6 < ?i-l , j *P i — 2 , j 
+ c 7<? i-l,j 2 + c 8“Pi , j 2 + c 9<Pi , j + c 10<Pi-l,j 
+ c ll<Pi,j+l + c 12 < Pi,j-l + c 13«Pi-2,j + c 12 r 

[15] R if j - civi, j + c 2 «Pi , j«Pi-l, j + 

+ c 4<Pi-l, jVi, j + c 5‘Pi,j < Pi-2,j + c 6 l Pi-l,j < Pi-2,j 
+ c 7 'P i — 1 , j 2 + c 8<Pi,j 2 + C 9<P i.j + c 10«Pi-l,j 
+ c 1 1*P i , j +1 + c 12 < Pi,j-l + c 13^i-2,j “ c ll r 
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Now, the previous 

expressions are differentiated with 

respect to <Pii ( jj to 

give the elements of the 

Jacob i an 

matrix 

< 3R i , j/ a< Pii , jj 

) This is achieved as follows 

For a 

general point (i,j). 


ii 

j j 

Derivative of 


i 

J-l 

c 1 2 


L 

J 

ci + c 3 <p i + 1> j + c 4 <p i 

-1 • j 

i 

J+l 

C 11 


i - 2 

j 

c 5 'P i + 1 , j + c 6Pi-l,j 

+ c 13 

i - 1 

j 

c 2*i+l,j + c 4«Pi,j + 

c 6Pi-2 , j 



+ 2c 7Vi _i,j + cio 


i + 1 

j 

c 2<Pi-l,j + c 3<Pi,j + 

c 5<Pi-2, j 



+ 2c 8 <Pi + i ( j + c 9 


ITE-1 

JB-2 

0 


ITE-1 

J B — 1 

0 


ITE-1 

JB 

0 


ITE-1 

JB + 1 

0 


ITE 

JB-2 

0 


ITE 

JB— 1 

0 


ITE 

JB 

0 


ITE 

JB + 1 

0 




For po ints [ 1 ] , 


9 1 


i i 

j j 

Derivative of Rj j 

i 

J -1 

c 1 2 

i 

J 

C 1 + c 3«Pi + l,j + c 4Vi-l,j + c xl 

l 

j+1 

0 

i - 2 

j 

c 5*i+l,j + c 6<Pi-l,j + c 1 3 

i- 1 

J 

c 2<Pi + l,j + c 4*i,j + c 6^i-2,j 

+ 2c 7 v i _ 1> j + c 10 

i + 1 

J 

c 2 < Pi-l,j + c 3^ i , j + c 5*i-2,j 
+ 2c 8 <p i + 1>J + eg 

ITE-1 

JB-2 

0 

ITE-1 

JB-1 

0 

ITE-1 

JB 

0 

ITE-1 

JB + 1 

0 

ITE 

JB-2 

0 

ITE 

JB-1 

0 

ITE 

JB 

0 

ITE 

JB + 1 

0 



For points [ 2 ] , 
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ii 

j j 

Derivative of R^ j 

i 

J -1 

0 

1 

J 

ci + c 3 <p i + 1 j + c 4 (p i _ liJ + c 12 

i 

j+1 

C 11 

i - 2 

j 

c 5^i + l,j + c 6<Pi-l,j + c 1 3 

i- 1 

j 

c 2^i + l , j + c 4 ! Pi , j + c 6<Pi-2,j 
+ 2 c 7 «pi_i ( j + c 10 

i + 1 

j 

c 2 i — 1 , j + ^i.j + c 5^i-2 , j 

+ 2c 8l p i + 1> j + eg 

ITE-1 

J B — 2 

0 

ITE-1 

J B — 1 

0 

ITE-1 

JB 

0 

ITE-1 

JB + 1 

0 

ITE 

J B — 2 

0 

ITE 

JB-1 

0 

ITE 

JB 

0 

ITE 

JB+1 

0 



For points [ 3 ] , 
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ii 

j j 

Derivative of j 

i 

J-l 

c 1 2 

i 

J 

ci + c 3 <p i + lf j + c u<Pi-l ,3 

i 

J + l 

C 11 

1-2 

j 

c 5^i + l , j + c 6 i - 1 , j + c 1 3 

i - 1 

j 

c 2'Pi + l , j + c 4'i p i , j + c 6 ( Pi-2 , j 
+ 2 c 7 <p i _ij + c 10 

i + 1 

j 

c 2Vi-l,l + c 3‘Pi,j + c 5 ^ i - 2 , j 
+ 2c 8! p i + 1( j + c 9 

ITE-1 

JB-2 

- CHC + 0.5T!) 

ITE-1 

J B — 1 

- cuC-1.51!) 

ITE-1 

JB 

- c 11 (+1.5T 1 ) 

ITE-1 

JB+1 

- cnC-o.sT].) 

ITE 

JB-2 

- C]_]_( + 0 5 T 2 ) 

ITE 

JB— 1 

- c n (-1.5T 2 ) 

ITE 

JB 

- C11 (+1.5T 2 ) 

ITE 

JB + 1 

- c x L ( -0 5T 2 ) 
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For points [ 4 ] , 


ii 

j j 

Derivative of 

i 

J-l 

c 1 2 

i 

j 

C 1 + c 3¥>i + l,j + c 4^i-l,j 

i 

J+l 

C 11 

i - 2 

j 

c 5^i + l , j + c 6*i-l,j + c 1 3 

i- 1 

J 

c 2^i + l , j + c 4<Pi,j + c 6'?i~2 , j 
+ 2c 7V ? i _ 1;J + c 10 

i + 1 

J 

c 2 < Pi-l,j + c 2Vi,j + c 5 i - 2 , j 
+ 2c 8 cp i + 1)J + c 9 

ITE-1 

JB-2 

+ cx2 ( +0 • 5T 1 ) 

ITE-1 

J B— 1 

+ c 12 (-1 . 5Tx) 

ITE-1 

JB 

+ c 12 ( + 1 • 5T 1 ) 

ITE-1 

JB + 1 

+ cj_2 ( -0 • 5T l) 

ITE 

JB-2 

+ c ]_ 2 ( + 0 5 T 2 ) 

ITE 

J B— 1 

+ C 12 (-1 • 5T2 ) 

ITE 

JB 

+ c^2 ( +1 • 5 t 2) 

ITE 

JB + 1 

+ c 2 ( “0 5 T 2 ) 
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For a subsonic free stream : 
For points [5], 


i i 

j j 

Derivative of 

i 

J -1 

c 1 2 

1 

J 

ci + c 4 <p i _ lj 

i 

J+l 

C 11 

i - 2 

J 

c 6<Pi-l,j + c 1 3 

i- 1 

J 

c 4<Pi j + c 6Vi _ 2>J 
+ c 10 

i + 1 

J 

0 

ITE-1 

JB-2 

0 

ITE-1 

J B — 1 

0 

ITE-1 

JB 

0 

ITE-1 

JB + 1 

0 

ITE 

JB-2 

0 

ITE 

J B— 1 

0 

ITE 

JB 

0 

ITE 

JB+1 

0 


2c 7^i- 1 , j 
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i i 

j j 

Derivative of Ri,j 

i 

J-l 

c 1 2 

i 

j 

ci + c 3 vi + l,j + c 4 < P i - 1 , j 

i 

J+l 

0 

i - 2 

J 

c 5Vi+l,j + c 6 ^i-l,j + c 1 3 

i - 1 

j 

c 2*i+l,j + c 4*i,j + c 6*i-2,j 

+ 2 c 7 ! Pi _ li j + cio 

l + l 

j 

c 2 < Pi-l,j + c 3<^i , j + c 5<Pi-2,j 
+ 2 cg<p i + i ( j + c 9 

ITE-1 

JB-2 

- c 11 (+0.5T 1 )/4 

ITE-1 

J B — 1 

- c 11 (-1.5T 1 )/4 

ITE-1 

JB 

- c 11 (+1.5T 1 )/4 

ITE-1 

JB + 1 

- c 11 (-0.5T 1 )/4 

ITE 

JB-2 

- cxx ( + 0 • 5 T 2 ) / 4 

ITE 

JB -1 

- C 11 (-1.5T 2 )/4 

ITE 

JB 

- C 11 (+1.5T 2 )/4 

ITE 

JB + 1 

- cxx(-0 . 5T 2 )/4 



For points [ 7 ] , 


ii 

j j 

Derivative of j 

i 

J-l 

c 1 2 

i 

j 

ci + c 3Vi + 1 j - c 4 r/2 

i 

J+l 

C 11 

i - 2 

J 

c 5*i+l,j ~ c 6 r/2 + c 13 

l-l 

j 

0 

i + 1 

J 

- c 2 r/2 + c 3 (p i( j + C 5<Pi_ 2>J 
+ 2c 8 <p i + 1> j + c 9 

ITE-1 

JB-2 

(* c 2V>i + l, j/2-c 4 <p ii j/2-c 6Vi _ 2 , j / 2 
+ c 7 r/2-c 10 /2) ( + 0 . 5 T L ) 

ITE-1 

J B — 1 

< -c 2*1+1 , j/2-c 4 ^i f j/2-c 6 Vi-2 , j/2 
+c 7 r/2-c 10 /2) (-1 5T X ) 

ITE-1 

JB 

< * c 2*i+l , j/2-c 4 <Pi f J /2-c 6 <p i _ 2 f j/2 

+c 7 r/2-c 10 /2) (+1 5T X ) 

ITE-1 

JB+1 

( -C2<Pi + l , j/2-c 4 <pi , j/2-c 6 <p i _ 2 f j/2 

+c 7 r/2-c 10 /2) (-0 5T L ) 

ITE 

JB-2 

( - c 2«Pi + l , j/2-c 4 <Pi f j /2-c 6 <ip 1 _2 , j/2 

+c 7 r/2-c 10 /2 ) ( +0 . 5T 2 ) 

ITE 

J B— 1 

(- c 2*i+l, j/2-c 4 <Pi f j/2-c 6l pi_2 t j/2 

+ c 7 r/2-c 10 /2) ( -1 . 5 T 2 ) 

ITE 

JB 

( -c 2 <Pi + l , j/ 2 - c 4<Pi , j/2-C6<Pi_ 2 ( j/2 
+c 7 r/2-c 10 /2) ( +1 . 5T 2 ) 

ITE 

JB + 1 

( - c 2*i+l , j/2-c 4 <Pi , j/2-c 6 <pi_ 2 , j/2 


+c 7 r/2-c 10 /2) (-0 . 5T 2 ) 
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For points 

[ 8 ] , 


i i 

jj 

Derivative of Ri j 

i 

J-l 

0 

i 

J 

ci + c 3 <p i + i ( j + C49> 1 -i tJ 

i 

J+l 

C I 1 

i - 2 

j 

c 5 <Pi + l,j + c 6 ( Pi-l , j + c 1 3 

i - 1 

J 

C 2<P i + l,j + c 4?i , j + c 6^i-2,j 



+ 2 c 7 cp i _i i j + cio 

i + 1 

j 

c 2<Pi-l,j + c 3 < Pi,j + c 5^i-2,j 



+ 2 c 8 Vi+l,j + c 9 

ITE-1 

JB-2 

- 3 c ! 2 (+0.5T X ) / 4 

ITE-1 

J B — 1 

- 3 c x 2 (-1.5T!) / 4 

ITE-1 

JB 

- 3 c x 2 ( + 1 . 5T]_ ) / 4 

ITE-1 

J B + l 

- 3 c x 2 (— 0 . 5T X ) / 4 

ITE 

JB-2 

- 3 c x 2 ( + 0 . 5T 2 ) / 4 

ITE 

JB-1 

- 3 c i 2 (-1 . 5T2 ) / 4 

ITE 

JB 

- 3 ci 2 (+1 • 5T2 ) / 4 

ITE 

J B + l 

- 3 c 12 ( — 0 . 5 T 2 ) / 4 



For points [ 9 J , 


ii 

i 

i 

i 

i - 2 
i- 1 

i + 1 
ITE-1 

ITE-1 

ITE-1 

ITE-1 

ITE 

ITE 

ITE 

ITE 


j j 

J-l 

j 

J+l 

j 


J 


j 

JB-2 


J B— 1 


JB 


JB + 1 


JB-2 


J B — 1 


JB 


JB + 1 


Derivative of j 
c 1 2 

C 1 “C3 r + c 4<Pi-l,j 
C 11 

- c 5 r + 6^ i— 1 , j + c 1 3 

- c 2 r + 4 ‘Pi , j + c 6‘Pi-2 , j 

+ 2c 7 ip£_ 1 t j + c 10 

0 

(-c 2 <Pi-i , j-C3<Pi , j-c 5 <p i _ 2 , j 

+ 2cgr— c 9 ) (+0 . 5T]_ ) 

(~c 2 <Pl-l , j-c 3 <Pi t j-C5<Pi_ 2 f j 
+2c 8 r-c 9 ) (-1 . 5T X ) 

(-c 2 9>i_i , j-C3<Pi f j-c 5 ^i_ 2 f j 

+ 2c 8 r-c 9 ) ( + 1 . 5T]_ ) 

<-C 2 «Pi_l t j-c 3 Vi , J-C 5 <Pi_ 2 _ J 

+ 2c 8 r-c 9 ) (-0 5 T x ) 

<-C2<Pi-l , j- c 3<Pi , j- c 5<Pi-2 , j 
+2c 8 T-c 9 ) (+0 . 5T 2 ) 

(-c 2 vi-l , j-c 3 «Pi ( j-c 5 <Pi_ 2 f j 

+2c 8 r-c 9 ) (-1 . 5T 2 ) 

(-c 2 Vi-l , j-c 3 SPi, j - 05 ^ 1-2 , j 
+2c 8 r-c 9 ) (+1.5T 2 ) 

<“ c 2^i-l , j“ c 3^i , j" c 5‘Pi-2 , j 
+ 2c 8 r-c 9 ) ( -0 . 5 T 2 ) 
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For po ints [ 10 ] , 


i i 

j j 

Derivative of 

R i 

. j 

i 

j-1 

c 1 2 




1 

J 

ci + 

c 4^i-l , j 



i 

J+l 

0 




i - 2 

J 

c 6*l- 

l.j + c 13 



i-1 

J 

c U < Pi , j + C 6 ( P±- 

2 , j 



+ c 10 




i + 1 

J 

0 




ITE-1 

JB-2 

- c 1X 

(+0. 5T X ) 

/ 

4 

ITE-1 

JB-1 

- c u 

(-1.5T!) 

/ 

4 

ITE-1 

JB 

- C11 

(+1 . 5T X ) 

/ 

4 

ITE-1 

JB + 1 

- cu 

( —0 . 5T X ) 

/ 

4 

ITE 

JB-2 

- c xl 

(+0 5T 2 ) 

/ 

4 

ITE 

JB-1 

“ C 11 

(-1. 5T 2 ) 

/ 

4 

ITE 

JB 

” C 11 

(+1 . 5T 2 ) 

/ 

4 

ITE 

JB+1 

" C 11 

( -0 . 5T 2 ) 

/ 

4 


2c 7<Pi-l , j 



- 

For points 

[ 11 ] . 

101 


i i 

j j 

Derivative of j 


i 

J-l 

c 12 


i 

J 

ci + C3V1+1.J ~ c 4 r / 2 

- 

i 

J+l 

0 


i - 2 

j 

C5<Pi + l,j “ C6T/2 + c 1 3 


1 - 1 

J 

0 


i + 1 

j 

- c 2 T + C 3 ^i j + c 5^1-2 , j 

- 

ITE-1 

JB-2 

+ 2 c 8 <Pi + l , j + c 9 

( -C 2 P 1 + 1 , j/2- c 4<Pi , j / 2 " c 6 ^ 1-2 , j /2 

- 

ITE-1 

J B— 1 

+ C7/2-cio/2-ch/4) ( + 0 5Ti ) 

( -C 2 <Pi + l , j/ 2 -c 4 <Pi f j/ 2 -c 6<Pi _2 , j /2 

- 

ITE-1 

JB 

+C 7 / 2 -C 10 / 2 -CH/ 4 ) (~1.5Ti) 
(-c 2 <Pi + l, j/ 2 - c 4 ( Pi, j/ 2 ~ c 6Vi-2 , j / 2 


ITE-1 

JB + 1 

+ C 7 / 2 -C 10 / 2 -C 11 / 4 ) ( + 1 5 T 1 ) 

( *c 2 <Pi + l > j/ 2 -c 4 < Pi ,j/ 2 _c 6^i-2 , j / 2 


ITE 

JB-2 

+ C7/2-cio/2-ch/4) (-0 5Ti) 

(-c 2 <Pi + l,j/ 2 - c 4 < Pi,j/ 2 _c 6 ( Pi-2,j / 2 

. 

ITE 

JB — 1 

+C7/2-cio/2-cii/4) (+0 5T 2 ) 

( -C2<Pi + l , j/2" c 4?i , j/ 2 “ c 6 ^i - 2 , j / 2 

- 

ITE 

JB 

+c 7 /2-cio/2-cii/4) (-1 • 5T 2 ) 

(-c 2 ‘Pi + l f j/2-C4«Pi ) j/2-C6(Pi_2,j/2 


ITE 

JB + 1 

+C7/2-cio/2-cii/4) (+1 . 5T 2 ) 

( -C2<Pi + l , j/2-C4<Pi , j / 2 ~ c 6^i-2 , j/ 2 




+C7/2-cio/2-ch/4) (-0.5T 2 ) 
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For 
i i 
1 
i 
i 

i - 2 
i-l 
1 + 1 

ITE 

ITE 

ITE 

ITE' 

ITE 

ITE 

ITE 

ITE 


points [ 12 ] , 

j j 

J-l 

J 

J+l 

j 

J 

j 


1 JB-2 


1 J B — 1 


1 JB 


1 JB + 1 


JB-2 


J B — 1 


JB 


JB+1 


Derivative of R^j 

0 

ci + C 3 P 1 + 1 J - c 4 r /2 
C 11 

cS'Pi + l.j “ c 6 r / 2 + c 1 3 
0 

- c 2 r /2 + c 3 cp i j j + c 5 <p i _ 2 ,j 

+ 2 cgtpi + i t j + c 9 

(~C 2 <Pi + l , j/2- c 4Pi , j/ 2 " c 6 ^i -2 , j / 2 
+C7r/2-cio/2-3ci2/^) (+0.5Ti) 
(-c 2 <Pi + l , j/ 2 -c 4 <Pi f j/2-c 6 cp i _ 2 , j/2 
+C7r/2-cio/2-3ci 2 / 4 ) (~ 1 - 5 T l) 
(-c 2 l Pi + l,j/ 2 -c 4 <Pi > j/ 2 -C 6 'Pi -2 , j / 2 
+ c7r/2-cio/2-3ci 2 /^) ( + 1.5T]_) 

(-e 2^1 + 1 , j/ 2 -c 4 ^i , j/2-c 6 ^i_ 2 , j /2 
+c 7 r/2-c 1 o/2-3c 12 /‘0 (-0 5T i ) 

(-c 2 Vi + i i j/ 2 -c 4 <Pi i j/ 2 -c 6 (p i _ 2 > j /2 

+C7r/2-cio/2-3ci 2 /^) ( + 0-5T 2 ) 
(“C 2 Vi + l , j/ 2 ~c 4 <Pi , j/2-c 6 <Pi-2 , j / 2 
+c 7 r/2-c 10 /2-3c 12 /^) (“1 • 5T 2 ) 
(-c 2 <Pi + l , j/2-c 4 <Pi ( j/2-c 6 ^i _ 2 , j/2 

+cjr/2-ciQ/2-3ci2/ 1 *) ( + 1.5T 2 ) 
(-c 2 <Pi + l, j/ 2 -c 4 9 i, j / 2 -C 6 Vi -2 , j / 2 
+C7r/2-cio/2-3ci 2 /^) (-0.5T 2 ) 



For points [13], 


i i 
i 
1 
i 

i - 2 
i - 1 

i + 1 
ITE-1 

ITE-1 

ITE-1 

ITE-1 

ITE 

ITE 

ITE 

ITE 


j 3 
J -1 
J 

J+l 

j 

j 

j 

JB-2 
J B — 1 
JB 

JB + 1 
JB-2 
J B — 1 
JB 

JB + 1 


Derivative of R.£ ( j 

0 

C 1 “ c 3 r + c 4‘Pi-l , j 
C 11 

- c 5 r + c 6 <Pi_i ( j + c 13 

- c 2 r + c a <Pi t j + c g(p i _2 , j 
+ 2c7<Pi-l,j + CIO 

0 

(-c 2 < Pi-l i j- c 3 < Pi,j -c 5 i Pi-2 ,j 
+ 2cgr-C9~3c]_ 2 /4) ( + 0 5T]_) 

(-c 2 <Pi_i t j-c 3 <Pi t j-c 5 <Pi_ 2 t j 
+2cgr-C9~3ci 2 /4) (-1 5Tj_) 

(-c 2 <Pi_i t j-c 3 <Pi f j-c 5 <Pi_ 2 , j 
+ 2cgr-C9-3c]_ 2 /4) ( + 1 . 5T]_ ) 

(-c 2 <Pi-i , j-c 3 <Pi , j- c 5^i-2 , j 
+ 2cgr-C9-3c]_ 2 /4) (-1.5T]_) 

(-c 2 <p i _i ( j-c 3 S £> iiJ -c 5 !p i _ 2i j 

+ 2 c gT-c 9 - 3 c ^ 2 /4 ) (+0.5T 2 ) 

(-c 2 Vi_i , j-c 3 <Pi ( j-c 5 <Pi_ 2 , j 
+2c 3 r-c9-3c^ 2 /4) (-1.5T 2 ) 
(-c 2 <Pi-l, j-c 3 «Pi , j"C5<Pi- 2 , j 
+2cgr-C9~3ci 2 /4) (+1.5T 2 ) 
(-c 2 <Pi-l , j"C3<Pi , j"C5^i — 2 , j 

+ 2 cgT — 09 — 3c^ 2 /4) ( —0 . 5 T 2 ) 
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For points [ 14 ] , 


i i 

j j 

i 

J -1 

i 

J 

i 

j+1 

1-2 

j 

i - 1 

j 

i + 1 

j 

ITE-1 

JB-2 

ITE-1 

JB-1 

ITE-1 

JB 

ITE-1 

JB + 1 

ITE 

JB-2 

ITE 

JB-1 

ITE 

JB 

ITE 

JB + 1 


Derivative of R< 

i 

c 1 2 

ci + c 4 ,p L _i i j 
C 11 

c 6 ( Pi-l , j + c 13 
c 4 < Pi , j + c 6^i-2,j 

+ c 10 
0 

c 12 (+0.5T1) 
c i 2 (-1.5T1) 
c i 2 (+1 5 Ti) 
ci2 ( -0 . 5 T ! ) 
c 12 ( + 0 • 5 t 2 ) 
ci2 (-1.5T2) 
c 12 ( + 1 5T2 ) 
ci2 (-0.5T2) 


+ 2c 7 cp i _ ljJ 



Fo r points [15], 


i i 
i 
i 
i 

i - 2 
i - 1 

i + 1 
ITE-1 

ITE-1 

ITE-1 

ITE-1 

ITE 

ITE 

ITE 

ITE 


j j 

J-l 

J 

J+l 

j 

j 

j 

JB-2 
J B — 1 
JB 

JB+1 

JB-2 

JB-1 

JB 

JB+1 


Derivative of Rf , j 
c 1 2 

c 1 - c 3 r + C 4 
C 11 

- c 5 r + c 6 ^i_ij + c 13 

- c 2 r + c 4‘P i , j + c 6<Pi-2,j 
+ 2 c 7 4 > i _ 1> j + c 10 

0 

(-c 2 Vi-l, j-c 3 vi, j-c 5 <Pi- 2 ,j 
+ 2c g T-c 9 -c 1 1 ) (+0.5T]_) 
(-c 2 ?>i-l , j-C3^i, j-c 5 Vi- 2 , j 

+2cgr-C9-cxi) (-1 . 5T^) 

(-c 2 Vi-l , j-c3<Pi t } ~c 5 <Pi-2 , j 
+ 2 cg T-c 9-c x i ) ( + 1 . 5 T]_ ) 

(~C 2 <Pi-l , j-C 3 <Pi , j"C5(Pi_ 2 ( j 

+ 2c8r-C9-cn ) (-0.5Tx) 

(-c 2 <Pi-l , j“C 3 ^i , j — ° 5<p i _ 2 , j 
+2csr-C9-cxx) (+0 . 5T 2 ) 
(-c 2 <p i _ 1 ( j-c 3<Pi t j-c 5 v i _ 2 ( j 
+ 2cgr-C9-cxi) ( -1 • 5 T 2 ) 
<" c 2 ?i-l , j ~ c 3 St?i , j- c 5^i-2 , j 
+ 2 c 3 T-c 9 -c x i ) ( + 1 . 5T 2 ) 
(-c 2<Pi _x , j-c 3 <Pi f j-C 5 <Pi _ 2 _ j 
+ 2 c g T-c 9 -c x x ) ( -0 . 5T 2 ) 
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For a supersonic free stream : 
For points [ 5 ] , 


i i 

j j 

Derivative of j 

i 

J-l 

c 1 2 

i 

J 

C 1 + c 2'Pi-l , j + 2c 3<Pi,j 
+ c 4‘Pi-l,j + c 5^i-2 , j + 
+ eg 

1 

J+l 

C 11 

i- 2 

j 

c 5 P i , j + c 6<P i — 1 , j + c 1 3 

i- 1 

J 

(c 2 + c4) ( p iiJ + c 6 'p x _ 2 , j 
+ 2c 7 <p i _ 1> j + c 10 

i + 1 

j 

0 

ITE-1 

JB-2 

0 

ITE-1 

J B— 1 

0 

ITE-1 

JB 

0 

ITE-1 

JB + 1 

0 

ITE 

JB-2 

0 

ITE 

J B— 1 

0 

ITE 

JB 

0 

ITE 

JB + 1 

0 


2c 8^i , j 
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For po ints 

[6] , 


i i 

j j 

Derivative of j 

i 

J-l 

c 1 2 

1 

j 

ci + c 3< p 1 + 1( j + c 4¥ > i _ 1 j 

i 

j+1 

0 

i - 2 

J 

c 5^i + l , j + c 6<Pi_i,j + C 13 

1-1 

J 

c 2*i+l,j + c 4 <p i(J + C6<Pi_ 2p j 
+ 2c 7 <p i _ li j + cio 

1 + 1 

j 

c 2<Pi-l,j + c 3<Pi,j + c 5'Pi-2,j 
+ 2c 8 cp i + 1( j + eg 

ITE-1 

JB-2 

0 

ITE-1 

J B — 1 

0 

ITE-1 

JB 

0 

ITE-1 

JB + 1 

0 

ITE 

JB-2 

0 

ITE 

J B— 1 

0 

ITE 

JB 

0 

ITE 

JB+1 

0 



For points 

[7] . 


i i 

j j 

Derivative of 

1 

J-l 

c 1 2 

i 

j 

ci + c 3 <p i + 1>J 

i 

J+l 

C 11 

i - 2 

j 

c 5 < Pi + l , j + c 1 3 

i - 1 

j 

0 

1 + 1 

j 

c 3^i , j + c 51Pi- 



+ c 9 

ITE-1 

JB-2 

0 

ITE-1 

JB-1 

0 

ITE-1 

JB 

0 

ITE-1 

JB + 1 

0 

ITE 

JB-2 

0 

ITE 

JB-1 

0 

ITE 

JB 

0 

ITE 

JB + 1 

0 


+ 2c 8! p i + li 
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For points 

[8] , 


i i 

j j 

Derivative of j 

i 

J-l 

0 

i 

J 

ci + c 3< p i + 1>J + C4*>i_l (J 

i 

j+1 

C 11 

i - 2 

J 

c 5^i+l,j + c 6 ‘P i — 1 , j + c 1 3 

i - 1 

J 

c 2*i+l,j + c 4'Pi,j + c 6<Pi-2,j 
+ 2c 7 cp i _ 1> j + cxo 

1 + 1 

j 

c 2'Pi-l,j + c 3^i , j + c 5 ^ i — 2 , j 
+ 2c 8 <p i + liJ + c 9 

ITE-1 

JB-2 

0 

ITE-1 

JB-1 

0 

ITE-1 

JB 

0 

ITE-1 

JB + 1 

0 

ITE 

JB-2 

0 

ITE 

JB-1 

0 

ITE 

JB 

0 

ITE 

JB+1 

0 



For points [9], 


i i 

j j 

Derivative of j 

i 

J -1 

c 1 2 

i 

j 

C 1 + c 2^i — 1 , j + 2c 3'Pi 1 j 



+ + c 5 cp i _2 i j + 



+ eg 

i 

J+l 

C 11 

CM 

i 

H 

J 

C 5<P i,j + c 6*i-l,j + ci3 

i - 1 

J 

(c 2 + c4) Vi f j + c 6 cpi_ 2 , j 



+ 2c 7 9i-l,j + c 10 

i + 1 

J 

0 

ITE-1 

JB-2 

0 

ITE-1 

JB-1 

0 

ITE-1 

JB 

0 

ITE-1 

JB + 1 

0 

ITE 

JB-2 

0 

ITE 

JB-1 

0 

ITE 

JB 

0 

ITE 

JB+1 

0 


2c8<P! 



Ill 


For points 

[10] , 


i i 

j j 

Derivative of Rj_ ( j 

i 

J -1 

c 1 2 

i 

J 

C 1 + c 2 ‘P i - 1 , j + 2 c 3<p 
+ c 4 <p i _ 1) j + c 5 <pi- 2 , 
+ eg 

i 

J+l 

0 

i - 2 

J 

^i + l.j + c 6'Pi-l,j 

i - 1 

j 

c 2 < Pi + l , J + c U c ?l , j + 
+ 2c 7 <p i _i i j + cio 

i + 1 

j 

0 

ITE-1 

JB-2 

0 

ITE-1 

JB-1 

0 

ITE-1 

JB 

0 

ITE-1 

JB + 1 

0 

ITE 

JB-2 

0 

ITE 

JB-1 

0 

ITE 

JB 

0 

ITE 

JB + 1 

0 


i . J 

j + 2 c 8 <pj_ i j 

+ C13 
c 6 ^ i- 2 , j 



For points 

[11] , 


i i 

j j 

Derivative of 

i 

J -1 

c 1 2 

1 

J 

ci + c2<Pi + l t3 

i 

J+l 

0 

i - 2 

J 

c 5 'P i + l , j + c 1 3 

i - 1 

J 

0 

i + 1 

J 

C3^i,j + c 5 , Pi- 



+ eg 

HE — 1 

JB-2 

0 

ITE-1 

JB-1 

0 

ITE-1 

JB 

0 

ITE-1 

JB + 1 

0 

ITE 

JB-2 

0 

ITE 

JB-1 

0 

ITE 

JB 

0 

ITE 

JB + 1 

0 


+ 2c 8 ^ i+ i, 



For points 

[12] , 


i i 

j j 

Derivative of 

i 

J -1 

0 

i 

J 

ci + C3<p i + i (J 

1 

J+l 

C 11 

i - 2 

J 

c 5 ( Pi + l , j + c 13 

i - 1 

J 

0 

i + 1 

J 

c 3 Vi , j + c 5<Pi- 
+ eg 

ITE-1 

JB-2 

0 

ITE-1 

JB-1 

0 

ITE-1 

JB 

0 

ITE-1 

JB + 1 

0 

ITE 

JB-2 

0 

ITE 

JB-1 

0 

ITE 

JB 

0 

ITE 

JB+1 

0 


+ 2c 8 <p i + li 



For points [ 13 ] , 


i i 

j j 

Derivative of R.£ j 

i 

J -1 

0 

i 

j 

C I + cavi-i.j + 2c 3 v ifJ 

+ c 4‘Pi-l,j + c 5^i-2,j + 2c 8 ( P 
+ c 9 

i 

J+l 

C 11 

i- 2 

j 

c 5'Pi + l , j + c 6<Pi-l,j + c 1 3 

i - 1 

j 

c 2<Pi + l,j + c 4^ i , j + c 6^i-2 , j 
+ 2c 7 ipi_]_ i j + cio 

i + 1 

j 

0 

ITE-1 

JB-2 

0 

ITE-1 

J B — 1 

0 

ITE-1 

JB 

0 

ITE-1 

JB + 1 

0 

ITE 

JB-2 

0 

ITE 

JB-1 

0 

ITE 

JB 

0 

ITE 

JB+1 

0 



For points [ 14 ] , 


i i 

j j 

Derivative of j 

i 

J -1 

c 1 2 

i 

J 

C 1 + c 2 ^i-l,j + 2 c 3 Vi,j 
+ c 4 <p i _ 1>J + C 5 <p i _ 2(J + 2 c g<p 
+ c 9 

i 

J+l 

C 11 

i -2 

J 

c 5^i + l,j + c 6 ( Pi-l , j + c 1 3 

i - 1 

J 

c 2*i+l,j + c 4<Pi,j + C6Vi-2,j 
+ 2 c 7 <p i _ 1) j + cio 

i + 1 

j 

0 

ITE-1 

JB -2 

c 12 ( + 0 . 5 T ! ) 

ITE-1 

JB-1 

c 12 (-1 5T x ) 

ITE-1 

JB 

c 12 ( + 1 . 5T X ) 

ITE-1 

JB + 1 

c 12 (-0.5T!) 

ITE 

JB-2 

<=12 (+ 0 . 5T2 ) 

ITE 

JB-1 

c 12 (-1.5T 2 ) 

ITE 

JB 

ci 2 ( +1 . 5T 2 ) 

ITE 

JB + 1 

c 12 5 t 2 ) 
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For points [ 15 ] , 

ii jj 

i J-l 

i- J 


Derivative of R x j 
c 12 

c 1 + c 2 <Pi-l,j + 2 c ivi '2 
+ C 4 <p i _]_j + c 5 <p i _ 2 iJ + 2 cg (p it] 
+ c 9 


i j +1 


C 11 


J c 5Pi+l,j + c 6 ^i-l,j + c 1 3 

j c 2‘i t> i + l , j + c 4 l Pi , j + 

+ 2 c 7 <p i _ 1 j + c 10 


i + 1 

j 

0 



ITE-1 

JB-2 

- 

C 11 

(+0. 5T X ) 

ITE-1 

JB-1 

- 

C 11 

(-1 • 5T X ) 

ITE-1 

JB 

- 

C 11 

(+1. 5T X ) 

ITE-1 

JB + 1 

- 

C 11 

(-0 . 5T X ) 

ITE 

JB-2 

- 

C 11 

(+0. 5T 2 ) 

ITE 

JB-1 

- 

C 11 

(-1.5T 2 ) 

ITE 

JB 

- 

C 11 

(+1 . 5T 2 ) 

ITE 

JB+1 

- 

C 11 

(-0. 5T 2 ) 



The derivates of the residual with respect to the 
design variables (i e right hand sides of Eq (1)) are 
given by , 


XD t 


Points 


Derivates 


T 


M« 


[1] (+cnAfj/gj + 4 ) F ' i 

[2] (-C]_2drj/ gj ) F' u 

General + c 2 ' <Pi + l , j Vi-1 , j 

+ c 3'^i + l, jVi.j + c 4'<Pi-l , j<Pi , j 
+ c 5 '*i + l , j^i-2 , j + c 6 >i-l , j<Pi-2 , j 
+ c 7 '<p i-l,j 2 + c 8 > i + li j 2 + c 9 '<p i + 1) 

+ c 10 '<Pi_l,j + c 1 3'9 i _2 (J 


a 

[ 1 ] 

(-cuAif/gj+j,) 



[ 2 ] 

(+c 12 Af?/g J _^) 


C 

[ 1 ] 

(+C 11 Ar ? /gj +h ) 

F'l 


[ 2 ] 

( -c 2 A»?/ gj — 4 ) 

F'u 

L 

[ 1 ] 

C+CiiArj/gj + 4 ) 

F'l 


[ 2 ] 

(-ci 2 A'7/g J -h) 

F ' 
r u 

Note 

that the primes 

(') denote the 

partial derivative 


with respect to the design variable (XDj_). 
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